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In a recent article, Griinebohm et al. [Phys. Rev. B 84 132105 (2011)] report that they fail to 
reproduce the A2 U ferroelectric instability of Ti0 2 in the rutile structure calculated with density 
functional theory within the PBE-GGA approximation by Montanari et al. [Chem. Phys. Lett 
364, 528 (2002)]. We demonstrate that this disagreement arises from an erroneous treatment of 
Ti 3s and 3p semi-core electrons as core in their calculations. Fortuitously the effect of the frozen 
semi-core pseudopotential cancels the phonon instability of the PBE exchange-correlation, and the 
combination yields phonon frequencies similar to the LDA harmonic values. 

Griinebohm et al. also attempted and failed to reproduce the soft acoustic phonon mode instability 
under (110) strain reported by Mitev et al. [Phys. Rev. B 81 134303 (2010)]. For this mode the 
combination of PBE-GGA and frozen semi-core yields a small imaginary frequency of 9.8i. The 
failure of Griinebohm et al. to find this mode probably arose from numerical limitations of the 
geometry optimization approach in the presence of a shallow double well potential; the optimization 
method is not suitable for locating such instabilities. 

PACS numbers: 63.20.D-, 71. 15. Mb, 77.84.-s, 77.22.-d, 77.80.bn 



I. INTRODUCTION 

Rutile Ti0 2 is an incipient ferroelectric material. Over 
a 20 year period density functional theory has been used 
to develop an understanding of its dielectric response and 
lattice dynamical properties. It has also established the 
exquisite sensitivity of these properties to strain, volume 
changes and thus to approximations such as the treat- 
ment of electronic exchange and correlation^. 

In 2002 two of us demonstrated 2 that within the 
Perdew Burke and Enzerhof generalized gradient approx- 
imation (PBE-GGA)^ the A 2tl mode is unstable at the 
equilibrium geometry. The resulting structural instabil- 
ity corresponds to a ferroelectric distortion of the crystal 
in qualitative disagreement with the observed properties 
of rutile TiO^ 4 -. In the local density approximation 
(LDA)£ the equilibrium cell volume is slightly smaller 
and a stable structure is obtained, giving a quantitatively 
correct description of the lattice dynamics and dielectric 
properties. This motivated the choice of LDA functional 
for our study of Ref. H on the influence of strain on the 
ferroelectric properties and phonon modes of rutile Ti0 2 . 
The major conclusion of that work was the existence of 
a low frequency acoustic phonon branch over a broad 
region of the Brillouin Zone around (1/2,1/2,1/4) which 
became soft under negative pressure or anisotropic strain. 

In article Ref. 0, Griinebohm, Ederer, and Entel, pub- 
lished another study on the same topic, in which they 
attempted and failed to reproduce those central findings 
of our work. Specifically they report that that the PBE- 



GGA does not predict a phonon instability of the rutile 
structure, that there is no softening of the A 2n mode un- 
der (110) strain, nor is there a destabilization under (110) 
strain at (1/2,1/2,1/4) by a soft acoustic mode. 

It would be highly unsatisfactory for this discrepancy 
to remain unresolved in the literature; it would weaken 
the basis for future studies on this topic, and under- 
mine confidence in the reproducibility of DFT plane- wave 
methods for this type of study. Griinebohm et al. offer a 
speculation that the discrepancy arises from "the differ- 
ent potentials used and the resulting difference in lattice 
parameters" . We tested that hypothesis, and confirm it 
as the most likely origin of the discrepancy. We show 
that the error lies in the choice of pseudopotentials used 
by Griinebohm et al. specifically their treatment of the 
Ti 3s and 3p semicore electrons as part of the frozen 
pseudopotential core instead of as valence. 

Textbook wisdom holds that in the case of titanate 
perovskites the frozen semicore approximation introduces 
significant errors in ferroelectric properties and that 3p 
and possibly 3s states must be explicitly included^. In- 
deed most pseudopotential studies of Ti0 2 and titanate 
perovskites published since the mid 1990s have treated 
the Ti 3s and 3p electrons as valence, claiming that the 
frozen semicore approximation is inaccurate. However 
there is no report containing quantitative evidence of the 
consequences for phonon frequencies. 

Ghosez et a L 13 ' 14 performed band-by band decompo- 
sition of the Born effective charges and showed that the 
Ti 3p orbitals contribute a value of -0.22e to the effective 
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TABLE I. Pseudopotential configurations used. A single ul- 
trasoft projector was used for each core or semi-core state. 
Two ultrasoft projectors were used for each valence state ex- 
cept for UFl where a single projector was used for Ti 4s and 
UF3 where a single norm-conserving projector was used for 
Ti 4s. For all frozen core and semi-core states a pseudized 
core charge^ was used. 



Label UFl 




UF2 


UF3 


Ti ref. 2s z 2p b 3s / 


3p b 4s^ 


3s^3p b 4s^ 


4s^3(f 
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2s 2 2p A 


2s 2 2p 4 


r c /a 0.8 
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1.3 


nnncr/ao 0.5 
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Cutoff 1300 
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TABLE II. Calculated structural parameters of rutile Ti0 2 




a/A 


c/A 


u 


PW-LDA (UFl) 


4.550 


2.919 


0.3039 


PW-LDA (UF2) 


4.549 


2.919 


0.3039 


PW-LDA (UF3) 


4.597 


2.902 


0.3023 


LCAO-LDA^ 


4.548 


2.944 


0.305 


FP-LAPW-LDA^ 


4.558 


2.920 


0.3039 


PW-GGA (UFl) 


4.642 


2.963 


0.3051 


PW-GGA (UF2) 


4.642 


2.964 


0.3051 


PW-GGA (UF3) 


4.684 


2.953 


0.3035 


PW-GGA^ 


4.641 


2.966 


0.305 


PW-GGA^ 


4.664 


2.969 


0.3047 


lcao-ggV 


4.623 


2.987 


0.306 











a Ref. and QJj 
b Ref. 
c Ref. 
d Ref. 
c Ref. 



charge on the Ti atoms in SrTi0 3 and -0.43e in the case 
of BaTi0 3 perovskites. The same authors also argued 
that small contributions to this effective charge are re- 
sponsible for the destabilization of the cubic phase by 
the softening of an optic phonon and a phase transition 
to the rhombohedral phase. Tegner et al. 1 ® performed 
a comprehensive investigation on the effect of inclusion 
of Ti 3s and 3p on cohesive energies, phase stability and 
electronic density of states of metallic Ti, concluding that 
frozen semicore Ti pseudopotentials result in poor trans- 
ferability and that the relative stability of the u and hep 
phases is reversed. Deskinsii found an 0.3 eV difference 
in the cohesive energy of Ti0 2 and again that semicore 
Ti pseudopotentials exhibit worse transferability in the 
case of changes of oxidation state, with errors of up to 
0.1 eV in dissociation reaction energies of TiCl ra and ab- 
sorption energetics of a variety of molecules and atoms 
on rutile Ti0 2 surfaces. Holzwarth et al*& highlighted 



a large error in the atomic configuration energy for s-d 
promotion in Ti atoms of 90 meV obtained within the 
frozen-semicore approximation. However Perron et al^ 
showed that the effect of a frozen-core pseudopotential 
on structural parameters of Ti0 2 rutile and anatase is 
less than 0.5% provided that a nonlinear core correction 
term is included 2 ^. 

The question of how many electrons to treat as va- 
lence is independent of other aspects of the pseudopo- 
tential construction or formalism, and holds equally for 
norm-conserving, ultrasoft or PAW pseudopotentials 18 , 
as it concerns the inclusion or omission of the physics of 
polarizability of the Ti 3s and 3p electrons. 

In summary, previous investigations are consistent 
with a picture in which structural parameters are little 
affected by the treatment of 3s and 3p, relative energies 
of phases and cohesive energies are noticeably affected. 
Frozen semicore pseudopotentials have poor transferabil- 
ity resulting in significant errors in dielectric properties 
and chemical reaction energetics. As no results have been 
reported on the consequence for phonon frequencies and 
mechanical stability, we present them here. 



II. GGA INSTABILITY OF RUTILE 
STRUCTURE 

We performed calculations on bulk rutile Ti0 2 using 
the CASTEP code^i, with a series of pseudopotentials to 
test the effect of the frozen-core approximation on succes- 
sive shells. The pseudopotentials were of the ultrasoft va- 
riety generated using Vanderbilt's method 2 - 2 *. Small core 
radii were chosen to generate accurate but "hard" pseu- 
dopotentials (i.e. which require a large cutoff) for both 
elements to leave the frozen-core approximation as the 
dominant contribution to the error. Details of these po- 
tentials are listed in table |U Three sets of pseudopoten- 
tials are designated UFl, UF2, and UF3 denoting that 
the Ti Is, 2s and 2p, and 3s and 3p shells were frozen re- 
spectively. The UF3 set corresponds to the "large core" 
set used by Grunebohm et al. and the UF2 set corre- 
sponds to the "small core" set used by Montanari and 
Harrison and by previous authors. The UFl set gives a 
"nearly all electron" calculation in which only the Ti Is 
orbitals are in the core, and Ti 2s and 2p and O Is are 
explicitly treated as valence. The computational cost of 
treating core states as valence is surprisingly reasonable 
because the sharply-peaked core orbitals are represented 
mostly by the pseudized augmentation functions of the 
Vanderbilt scheme. A very fine FFT grid is required 
only for the electron density to represent the augmenta- 
tion charge, not for the soft orbital functions. This "all 
electron pseudopotential" approach has previously used 
successfully for calculations on iron at TPa pressures 2 ^. 

In each case a geometry optimization was performed 
to optimize both lattice parameters and the internal co- 
ordinate, followed by a finite displacement calculation 
of the dynamical matrix. This was diagonalized to yield 
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TABLE III. Calculated and experimental F-point phonon frequencies of bulk rutile Ti0 2 (cm 1 ). (TO frequencies only.) 



Mode LDA-UF1 LDA-UF2 LDA^ LDA-LAPW^ PBE-UF1 PBE-UF2 PBE^ PBE-LW PBE-UF3 Expt^ 



B2« 


99.2 


101.7 


104.0 


111.5 


31.0z 


28Ai 


79.2 


96.1 


91.2 


113 


A 2il 


125.1 


128.7 


154.4 


153.0 


89.6i 


106.6i 


86.3i 


143.7 


115.7 


142 




135.2 


138.8 


191.4 


156.4 


69.8i 


62.5i 


124.0 


137.7 


111.2 


189 




136.9 


137.1 


137.0 


138.2 


127.1 


150.5 


154.2 


131.2 


144.4 


142 


E u 


381.5 


382.4 


383.9 


384.8 


353.8 


355.1 


353.5 


383.2 


366.4 


388 




390.9 


392.4 


393.0 


398.2 


362.3 


356.4 


357.5 


400.5 


380.3 


406 




417.6 


419.5 


421.7 


416.0 


408.5 


417.9 


423.6 


407.6 


405.6 




E 9 


463.1 


463.5 


463.2 


463.7 


429.1 


428.9 


429.2 


472.2 


464.3 


455 


E« 


486.7 


487.3 


488.4 


487.0 


469.6 


469.7 


488.4 


480.0 


491.4 


494 


Alg 


609.9 


610.4 


611.6 


612.0 


567.8 


567.3 


565.9 


626.9 


611.1 


610 


Big 


817.3 


818.2 


824.7 


816.1 


770.8 


771.4 


774.3 


815.2 


794.6 


825 



a Ref.H 

b Ref.Qj]. These are pseudopotential-LAPW, not all-electron 

FP-LAPW calculations 
c Calculated using the optimized lattice parameters of UF2 



the gamma-point phonon frequencies. The Brillouin zone 
was sampled using a Monkhorst Pack 24 grid of dimen- 
sions 4x4x6 and plane-wave cutoffs for each USP set 
are listed in table Q] These choices yield a convergence of 
better than 1 meV/atom in total energy, 0.00025 A in lat- 
tice parameter. The geometry was relaxed until the force 
residual on each atom was less than 0.005 eV/A. For most 
phonons the frequencies were converged to better than 
0.1 cm , for the most sensitive phonon modes (which 
become imaginary) the maximum error was 3 cm -1 . 

The calculated structural parameters are shown in ta- 
ble |nl LDA lattice parameters for UF1 and UF2 are 
within 0.2% of those calculated using the reference all- 
electron techniques'^. Freezing the Ti semicore states 
slightly increases the lattice parameters but by less than 
1%, consistent with the values reported by Grimebohm 
et al. and by previous studie s 19 ' 20 . 

Phonon frequencies are presented in table IIII1 This 
shows clearly that (a) the rutile structure is predicted 
to suffer from a phonon instability within PBE and (b) 
that freezing the semicore states has the effect of elimi- 
nating the instability and yields all real frequencies. As 
the atomic reference configuration is the same as used by 
Grimebohm et al. it is highly likely that this explains the 
failure to reproduce the instability in their calculation. 
The change in frequencies when the UF3 pseudopoten- 
tials are used at the lattice parameters of the UF2 calcu- 
lation is much smaller than the UF2-UF3 change. This 
proves that the discrepancy between Griinebohm's result 
and ours is a direct consequence of the omission of elec- 
tronic polarizability and hybridisation of the Ti 3s and 
Sp states and not an indirect consequence of the volume 
change. 

The accuracy of the pseudopotentials is confirmed by 
the LDA results which are in excellent agreement with 
previous work including FP-LAPW calculations. The 
similarity of the UF2 results to the nearly all-electron 
UF1 values for both LDA and PBE demonstrates that the 



TABLE IV. Frequencies of 


the lowest 


mode at 


q=(0.5,0.5,0.25) calculated using 


using a y2 


x x 4 


supercell (cm - ). 






UF1 


UF2 


UF3 


PBE 77.1i 


72.6i 


9.8i 


LDA 50.0 


55.1 


74.8 



UF2 configuration contains all of the electronic freedom 
required for accurate lattice dynamics. This confirms the 
correctness of the "small-core" pseudopotential results of 
Ref. [3 where it was demonstrated that all-electron LCAO 
methods gave the same A2 U instability as the plane- wave 
pseudopotential calculations. We note one discrepancy 
with the prior results of Ref. 0, namely that the B 2u and 
E n modes also exhibit an instability in addition to the 
A 2u instability previously described. 



III. MODE SOFTENING AT Q=(l/2,l/2,l/4) 

Turning to the soft acoustic mode frequencies, we per- 
formed supercell phonon calculations using a \/2x\/2 x 4 
supercell and diagonalized at the commensurate wavevec- 
tor q=(0.5,0.5,0.25). Frequencies of the lowest acoustic 
modes are presented in table IIVI Clearly the UF3 pseu- 
dopotential makes a substantial error in this mode, which 
is stiffened when the effect of semicore polarization is ex- 
cluded. As with the T -point phonons, there is little dif- 
ference between the unfrozen core UF2 compared with 
the nearly all-electron UF1 frequencies. The UF1 and 
UF2 LDA results are in good agreement with our earlier 
calculations in Ref. @. However our current calculations 
predict that within the PBE approximation this acoustic 
mode has imaginary frequency for all three pseudopoten- 
tial configurations considered. As with the Gamma point 
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results, we would expect the UF3 case to most closely re- 
produce the results of Ref. @- 

Gruncbohm et al. state that "we could not repro- 
duce the destabilization of the system along an acous- 
tic phonon mode under (110) strain found in Ref. 82£, 
for which an energy gain of several meV has been pre- 
dicted, even if we use a 2.y2 x 2.\/2 x 4 supercell, which 
is commensurable with the displacement pattern of this 
mode" . 

However the geometry optimization method as em- 
ployed by Griinebohm et al. is unsuitable for finding 
the soft displacement pattern. Our experience is that 
quasi- Newton type optimization methods usually fail any 
attempt to identify a shallow soft mode, or at best con- 
verge extremely slowly. Discovery of a single set of 
atomic displacements which slightly lowers the energy 
in a high-dimensional search space of otherwise stiff di- 
rections is a highly computationally challenging proce- 
dure. Several factors contribute to the difficulty. First is 
the high dimensionality of the search space, 576 for the 
\/2 x \/2 x 4 cell employed (which is actually 4 times larger 
than strictly necessary to represent this soft eigenvector) . 
Within this space the desired eigenvector is represented 
by a single direction. 

The second factor is the shape and depth of the anhar- 
monic potential, which has a well depth of only 2.5 meV 
(see Fig. 5 of Ref. H) and a tiny gradient along the mode. 
For small displacement this yields a typical force on each 
Ti of 0.01 eV/A, with a maximum of 0.03 eV/A. Conse- 
quently the configuration at the initial points along the 
optimization would be deemed already converged accord- 
ing to the tolerance of 0.01 eV/A used by Griinebohm et 
al. In our experience a tolerance of 0.001 eV/A or even 
smaller is essential in such circumstances. 

The third factor is the difficulty of finding a suitable 
starting configuration for the optimization. In a super- 
cell of perfect rutile the gradient in every direction is 
zero, so an initial perturbation is required to break the 
supercell symmetry. A randomized displacement pattern 
must also excite the stiff modes and yield a much larger 
increase in the energy than the decrease being sought. 
The resulting geometry optimization almost always fails 
to find the soft mode displacement. We tested this for 
the soft acoustic mode in rutile Ti0 2 , using the UF2 
pseudopotential, a force tolerance of 0.001 eV/A and as 
expected the optimization returned to the symmetric su- 
percell configuration. 

From the UF3 phonon frequency results in table ITVl we 
would expect an even shallower potential and even tinier 
gradient as a result of the additional stiffening resulting 
from the neglect of semi-core polarization. It is therefore 
doubly unsurprising that Griinebohm et al.'s attempt to 



search for the soft displacement using this method was 
unsuccessful. This is absence of evidence for the soft 
mode, not evidence of absence. 

In our experience the only reliable way to explore a soft 
mode using geometry optimization is to first perform a 
phonon calculation, and to identify the soft eigenvector 
by diagonalization of a dynamical matrix. This eigen- 
vector can be used as an initial perturbation from the 
symmetric supercell. The geometry optimization must 
use very tight force and energy tolerances, but is usually 
able to proceed without interference from the stiff modes. 

IV. CONCLUSION 

We have demonstrated that the finding of two of us 
that rutile Ti0 2 is mechanically unstable to a ferroelec- 
tric A 2ll soft mode in the PBE-GGA approximation is 
robust. Griinebohm et aL's contrary finding can be at- 
tributed to an error in the PAW pseudopotentials they 
used, namely the frozen-core approximation applied to 
the Ti 3s and 3p electrons. Our convergent series of 
pseudopotentials to a highly accurate "nearly all elec- 
tron" limit demonstrates that the effect of freezing Ti 3s 
and 3p is to stabilise the PBE rutile structure of Ti0 2 . 
In fact this error almost cancels the PBE instability to 
yield T point frequencies very similar to the LDA result 
but the cancellation is less complete elsewhere in the Bril- 
loun Zone. 

We further argue that the geometry optimisa- 
tion methods with convergence tolerances reported by 
Griinebohm et al. are not capable of identifying shallow 
soft modes in the Brillouin zone. Our result of a soft 
mode at q=(l/2,l/2,l/4) is also robustly confirmed. 

Thus the situation is that in the incipient ferroelectric 
rutile-Ti0 2 the LDA and PBE-GGA approximations to 
density functional theory yield qualitatively distinct de- 
scriptions. Within PBE-GGA, rutile Ti0 2 is predicted 
to be ferroelectric with an equilibrium geometry that is 
unstable with respect to displacement along a number of 
phonon eigenvectors^. Within the LDA, rutile- Ti0 2 is 
stable and a reasonable description of its dielectric re- 
sponse and lattice dynamics is obtained. 
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